install.packages(knitr)
install.packages(sf)
install.packages(mapsf)
install.packages(packcircles)
install.packages(cartogram)
install.packages(recmap)
install.packages(dplyr)
install.packages(cartogramR)
# install.packages("https://cran.r-project.org/src/contrib/Archive/cartogramR/cartogramR_1.0-1.tar.gz", repos = NULL, type = "source")

library(sf)
communes <-
  st_read(
    "https://raw.githubusercontent.com/transcarto/rcartograms/main/data/isere.geojson",
    quiet = TRUE
  ) %>% st_transform(2154)
data <-
  read.csv(
    "https://raw.githubusercontent.com/transcarto/rcartograms/main/data/popisrere.csv",
    dec = ","
  )
communes = merge(x = communes[, c("id", "name", "geometry")],
                 y = data[, c("id",
                              "pop2018",
                              "agri",
                              "art",
                              "cadr",
                              "interm",
                              "emp",
                              "ouvr",
                              "retr")],
                 by = "id")
isere = st_union(communes)

knitr::kable(communes[c(0:10),], row.names = F, digits = 1)

dots = communes
st_geometry(dots) <- st_centroid(sf::st_geometry(dots),of_largest_polygon = TRUE)
dots <- data.frame(dots$id, dots["pop2018"], st_coordinates(dots))
dots = dots[,c("dots.id","X","Y","pop2018")]
colnames(dots) <- c("id","x","y","v")
dots <- dots[!is.na(dots$v),]
knitr::kable(dots[c(0:5),], row.names = F, digits = 1)

library("packcircles")

k = 500 # pour ajuster la taille des cercles
itermax = 10 # nombre d'iterations

dat.init <- dots[,c("x","y","v")]
dat.init$v <- sqrt(dat.init$v * k)
simulation <- circleRepelLayout(x = dat.init, xysizecols = 1:3,
                                wrap = FALSE, sizetype = "radius",
                                maxiter = itermax, weights =1)$layout
knitr::kable(simulation[c(0:5),], row.names = F, digits = 1)

circles <- st_buffer(sf::st_as_sf(simulation, coords =c('x', 'y'),
                      crs = sf::st_crs(communes)), dist = simulation$radius)
circles$v = dots$v

template("Dorling (packcircles)", "maps/dorling1.png", scale = FALSE)
mf_map(isere, col = "#CCCCCC", border = "white", lwd = 0.5, add = TRUE)
mf_map(circles,col = col, border = "#6b4266", lwd = 0.5, add = TRUE)
dev.off()

library(cartogram)

library(cartogramR)

precarto <- precartogramR(communes, method = "GastnerSeguyMore")
summary(precarto)

GastnerSeguy <- cartogramR(communes, count="pop2018", method="GastnerSeguyMore", options=list(L=256, grid=TRUE, maxit = 5))

grid <- make_layer(GastnerSeguy, type = c("final_graticule"))

template("Gastner, Seguy & More", "maps/gastnerseguy.png", basemap = FALSE, scale = FALSE)
mf_map(GastnerSeguy$cartogram, col = col, border = "white", lwd = 1, add = TRUE)
mf_map(grid, col = NA, border = "#6b4266", lwd = 0.05, add = TRUE)
dev.off()

table = cbind(GastnerSeguy$initial_data[,c("id","pop2018")] %>% st_drop_geometry(),
      orig_area = GastnerSeguy$orig_area,
      final_area = GastnerSeguy$final_area,
      errors = residuals(GastnerSeguy, type = "relative error")*100
      )
knitr::kable(table[c(0:10),], row.names = F, digits = 1)

st_write(as.sf(GastnerSeguy),"data/GastnerSeguy.geojson", delete_layer=TRUE)

dotcartogram = function(x,var,itermax,onedot,radius){
crs = sf::st_crs(x)
coords <- st_coordinates(st_centroid(sf::st_geometry(x),of_largest_polygon = TRUE))
x <- x[c("id",var)] %>% st_drop_geometry()
x <- data.frame(x, coords)
colnames(x) <- c("id","v","x","y")
x$v <- round(x$v/onedot,0)
x <- x[x$v > 0,]
dots <- x[x$v == 1,c("x","y","v")]
rest <-  x[x$v  > 1,c("x","y","v")]

nb <- nrow(rest)
  for (i in 1:nb){
    new <- rest[i,]
    new$v <- 1
    for (j in 1:rest$v[i]){ dots <- rbind(dots,new)}
  }
  
dots$x <- jitter(dots$x)
dots$y <- jitter(dots$y)
dots$v <- radius
  
simulation <- circleRepelLayout(x = dots, xysizecols = 1:3,
                                  wrap = FALSE, sizetype = "radius",
                                  maxiter = itermax, weights =1)$layout

circles <- st_buffer(sf::st_as_sf(simulation, coords =c('x', 'y'),
                                    crs = crs), dist = radius) 
return(circles)
}

onedot = 1000
dc = dotcartogram(x = communes, var = "pop2018", itermax = 120,
                  onedot = onedot, radius = 600)

template("Dots Cartogram", "maps/dotcartogram.png", note = paste0("Un point =\n",onedot," personnes"), scale = FALSE)
mf_map(isere, col = "#CCCCCC", border = "white", lwd = 0.5, add = TRUE)
mf_map(dc, col = col, border = "#6b4266", lwd = 0.8, add = TRUE)
dev.off()

cartogram = st_read("https://raw.githubusercontent.com/transcarto/rcartograms/main/data/GastnerSeguy.geojson")
# Ou en local
# cartogram = st_read("data/GastnerSeguy.geojson")

grid = st_make_grid(cartogram, cellsize = 2000, square = FALSE)
grid = st_sf(id = rep("",length(grid)), geometry = grid)
ctr = st_centroid(grid)
ids = cartogram[,"id"] %>% st_drop_geometry()

for(i in 1:nrow(grid)){
  x = st_within(x = ctr[i,], y = cartogram, sparse = TRUE, prepared = TRUE)
  id = ids[unlist(x),][1]
  if (length(id) == 0){id = NA}
  grid[i,"id"] = id
}
grid <- grid[!is.na(grid$id),]

gridcartogram <- aggregate(x = grid, 
               by = list(grid$id),
               FUN = min)

varmax = sum(cartogram$pop2018)
nbcell = nrow(grid)
valcell = round(varmax / nbcell)

template("Template cartographique", "maps/hexcartogram.png", note = paste0("Un hexagone ≈\n",valcell," habitants"),basemap = FALSE, scale = FALSE)
mf_shadow(x = grid, col = "grey50", cex = 1, add = TRUE)
mf_map(grid, col = col, border = "white", lwd = 0.5, add = TRUE)
mf_map(gridcartogram, col = NA, border = "#6b4266", lwd = 1, add = TRUE)
dev.off()

c <- st_centroid(communes)
v <- st_voronoi(x = st_union(c))
v <- st_intersection(st_cast(v), st_union(communes))
v <- st_join(x = st_sf(v), y = c, join=st_intersects)
v <- st_cast(v, "MULTIPOLYGON")

template("Polygones de voronoi", "maps/voronoi.png", basemap = FALSE, scale = FALSE)
mf_map(x = v, col = col, border = "#6b4266", lwd = 0.5, add = TRUE)
dev.off()

VDougenik <- cartogram_cont(v, "pop2018", prepare = "none", itermax = 15, maxSizeError = 1.15)

template("Dougenik & voronoi", "maps/dvoronoi.png", basemap = FALSE, scale = FALSE)
mf_map(x = VDougenik, col = col, border = "#6b4266", lwd = 0.5, add = TRUE)
dev.off()

library(recmap)

sfc_as_cols <- function(x, geometry, names = c("x","y")) {
  if (missing(geometry)) {
    geometry <- sf::st_geometry(x)
  } else {
    geometry <- rlang::eval_tidy(enquo(geometry), x)
  }
  stopifnot(inherits(x,"sf") && inherits(geometry,"sfc_POINT"))
  ret <- sf::st_coordinates(geometry)
  ret <- tibble::as_tibble(ret)
  stopifnot(length(names) == ncol(ret))
  x <- x[ , !names(x) %in% names]
  ret <- setNames(ret,names)
  dplyr::bind_cols(x,ret)
}

coords <- st_centroid(communes)
coords <- st_transform(coords, crs="+proj=longlat +datum=WGS84 +ellps=WGS84")

df_recmap <- sfc_as_cols(coords)

rec_cartogram <- data.frame (x= df_recmap$x,
                             y=df_recmap$y,
                             # make the rectangles overlapping by correcting lines of longitude distance
                             dx = sqrt(df_recmap$pop2018) / 90 / abs((0.65 * 60 * cos(df_recmap$y*pi/180))),
                             dy = sqrt(df_recmap$pop2018) / 90 / (0.65 * 60),
                             z = sqrt(df_recmap$pop2018),
                             name = df_recmap$id)

template("RecMap", "maps/recmap1.png", basemap = FALSE, scale = FALSE)
plot.recmap(rec_cartogram, col = NA, border = col, lwd=4,  col.text = col)
dev.off()

# cartog <- recmap(df)
# template("RecMap", "maps/recmap2.png", basemap = FALSE, scale = FALSE)
# plot(cartog[!cartog$name %in% c("IS","MT"),],  col = col, border = "white")
# dev.off()

europe <-
  st_read(
    "https://raw.githubusercontent.com/transcarto/rcartograms/main/data/europe.geojson"
  ) %>% st_transform(3035)

coords = data.frame(st_coordinates(st_centroid(st_geometry(europe))))
bb <- lapply(st_geometry(europe), function(x) {
  st_bbox(x)
})
dx <- unlist(lapply(bb, function(x) {
  x[3] - x[1]
})) / 2
dy <- unlist(lapply(bb, function(x) {
  x[4] - x[2]
})) / 2

df <- data.frame(
  x = coords$X,
  y = coords$Y,
  dx = dx,
  dy = dy,
  z = europe$pop2008,
  name = europe$id
)

template("", "maps/recmap2.png", basemap = FALSE, scale = FALSE)
plot.recmap(
  df,
  col = NA,
  border = col,
  lwd = 4,
  col.text = col
)
dev.off()

cartog <- recmap(df)
template("", "maps/recmap3.png", basemap = FALSE, scale = FALSE)
plot(cartog[!cartog$name %in% c("IS", "MT"),],  col = col, border = "white")
dev.off()

communes$inverse = 1/communes$pop2018